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We study the dynamics of a granular gas heated by the stochastic thermostat. From a Boltzmann 
description, we derive the hydrodynamic equations for small perturbations around the stationary 
state that is reached in the long time limit. Transport coefficients are identified as Green-Kubo 
formulas obtaining explicit expressions as a function of the inelasticity and the spacial dimension. 



I. INTRODUCTION 

Granular assemblies -made up of macroscopic solid bodies undergoing dissipative interactions- have been shown 
P, [1[ to frequently exhibit flows similar to those of normal fluids, and for practical purposes are often described by 
phenomenological hydrodynamic equations, i.e. equations for the density, flow velocity and energy density. This is so 
even if energy is not a conserved variable (the condition is that the energy mode be a slow variable compared to the 
rest of excitations). When the dynamics of the grains can be partitioned into sequences of two-body collisions, there 
is support both from experiments and computer simulations for the usefulness of a Kinetic Theory description. This 
occurs, schematically, when there are no clusters nor jamming and the system remains fluidized Q. In the low density 
limit, the relevant dynamics is encoded in the one-particle distribution function, which obeys the inelastic Boltzmann 
equation [H, H| . This is the starting point for many of the formal derivations of hydrodynamic equations by applying 
similar tools and ideas as those used in the context of ordinary fluids @. In the free-cooling case the study of the 
existence and applicability of a hydrodynamic regime is rather complete for the inelastic hard sphere (IHS) model. 
The Navier-Stokes equations have been derived by the Chapman-Enskog expansion Q and also via the linearized 
Boltzmann equation [81410} , yielding equivalent Green-Kubo formulas for the transport coefficients . The problem 
for arbitrary densities has also been tackled in [12, [HJ applying linear response methods. Although the successful of 
the Navier-Stokes equations is remarkable, granular systems often require to go beyond this level of description. For 
these cases and close to a stationary state, a modification of the Chapman-Enskog expansion has been carried out 
taking into account rheological effects [3, EH . 

Some modifications of the IHS model have been introduced for the study of other situations, mainly when a 
stationary state is reached. For example, when in a vibrated system the stationar y st ate is quasi-homogeneous, or 
when the grains are immersed in an interstitial medium that acts as a thermostat [16H19I ] , the system can be effectively 
modelled as driven by some random energy source. The energy injection can be performed by applying a random 
force to each particle. Depending on the stochastic properties of this force, different kinds of thermostats are obtained 
[20| . one of the most used being the so-called stochastic thermostat, which consists of a white noise force acting on 
each grain [2ll - [34j |. This model has been less studied that its unforced version (free-cooling case). To our knowledge, 
the only derivations of the hydrodynamic equations are the one made in [2tjj which consider some variants of the 
present stochastic thermostat in which the heating may depend on the local temperature. The objective of this work 
is to go further in this direction and to derive the hydrodynamic equations for the actual homogeneous stochastic 
thermostat model in the low density limit. It is arguably the most commonly employed model in the simulations. At 
the Boltzmann equation level, we will consider states that are close to the homogeneous stationary regimes, thereby 
obtaining linear equations for the deviations of the hydrodynamic fields around their homogeneous counterparts, with 
explicit Green-Kubo formulas for the transport coefficients. Let us note that, very recently, a study in the same lines 



has been done in 3J| by a different method, applying the Chapmann-Enskog expansion to the Enskog equation. 

The plan of the paper is as follows. The model is first defined in section [TXJ We summarize the main properties 
of the most general hydrodynamic state through which the stationary regime is reached in the long time limit. This 
state was analyzed in detail in [35j and, as will be shown, plays an essential role in the hydrodynamic description 
for the present model. In section IIII1 the linearized Boltzmann equation is written and the relevant modes for the 
hydrodynamic description arc identified. These properties arc subsequently exploited in section IIVI to derive the 
linearized hydrodynamic equations. Finally, in section [Vj we present a short summary of the results obtained in the 
paper while details of the calculations are given in seven appendixes, at the end of the text. 
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II. THE MODEL 



Let us consider a dilute gas of N smooth inelastic hard spheres (d = 3) or disks (d = 2) of mass m and diameter 
a. These bodies collide inelastically with a coefficient of normal restitution a, independent of the relative velocity. 
If at time t there is a binary encounter between particles i and j, with velocities Vj(t) and V,-(i) respectively, the 
postcollisional velocities V£(t) and ~V'j(t) are 



V'i = V, - 

V' = V, + 



2 

1 + a 



(<7 • Vij)&, 

(<x ■ Vij)&, 



(1) 



where = V; — V,- is the relative velocity and <x is the unit vector pointing from the center of particle j to the 
center of particle i at contact. Between collisions, the system is heated uniformly by adding a random velocity to the 
velocity of each particle independently with certain frequency and with a given probability distribution. Let us define 
the jump distribution, P& t (Av), as the probability that a particle experiences a jump Av in the time interval At, 
that will be assumed to fulfill 



A lim J dy yj P At (y) = 0, lim ^ J ^P At (y)=e„ 2 , 



(2) 



where £q is the strength of the noise. If the variance of this distribution is small compared to the velocity scale in 
which the one-particle distribution, /(r, v, i), varies, the evolution equation is, in the low density limit, the Boltzmann- 
Fokkcr-Planck equation [22|, |36| 



d_ 

dt 



Vi 



d 



f(xi,t) = a 



d-l 



dx 2 S(r 12 )f (v l ,y 2 )f(x 1 ,t)f(x 2 ,t) + y^/^i,*). 



where we have introduced the field variable x = {r, v} and the binary collision operator To 

To(vi,v 2 ) = / do-e(v 12 • <x)(v 12 • <r)(a _2 6 o : 1 - 1). 



(3) 



(4) 



Here the operator b a 1 replaces the velocities vi and v 2 by the precollisional ones and v 2 given by 

1 + a 



Vl 

v 2 



2a 
1 + a 

2a 



(cr ■ v 12 )&. 



(5) 



Under the conditions noted above, the evolution equation does not depend on the details of the distribution Pa*, but 
only on its second moment, through the coefficient £q- 

It is convenient to introduce the hydrodynamic fields in the standard kinetic theory fashion, as the first velocity 
moments of the one-particle distribution function 

(6) 
(7) 

(8) 



n(r,t) = / dvf(r,v,t), 
n(r, t)u(r, t) = J dvvf(r, v, t), 
~n(r,t)T(r,t) = J dv^V 2 f(r, v, t), 



where V = v — u is the velocity of the particle relative to the local velocity flow, u. Let us also introduce the local 
thermal velocity through 



v(r,t) 



2T(r,i) 



1/2 



(9) 



We will see that, in terms of the homogeneous hydrodynamic fields, we can specify some relevant states at the 
Boltzmann equation level. It is known numerically that for a wide class of initial conditions, the system reaches 
a homogeneous stationary state [23[. Assuming that total momentum is zero, i.e. J dvv/(v,0) = 0, the state is 
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characterized by an isotropic stationary distribution, f s (v), which was studied in detail in 22[. There, the distribution 
was written as 



v 



(10) 



where n is the total density and v s is the stationary thermal velocity. We have also introduced the scaled distribution 
function, Xs, which is independent of the strength of the noise £ 2 , i.e. all the dependence of the distribution on £q 
is written in terms of the temperature. As this distribution is quite close to a Maxwellian, an expansion in terms of 
Sonine polynomials does make sense. In the so-called first Sonine approximation the function reads [22j 



Xs(c)^Xm(c)[1 + 4S 2 (c 2 )}, 
where xm is the Maxwellian distribution with unit temperature, S 2 (c 2 ) 



d(d+2) 
8 



(11) 

^4r-c 2 + ic 4 is the second So- 



nine polynomial, and a 2 is the kurtosis of the distribution. Within this approximation and neglecting non linear 
contributions in a 2 , the distribution function can be calculated as (22j 

s(a) = 16(l-a)(l-2a 2 ) 

2[ ' 73 + 56d - 24da - 105a + 30(1 - a)a 2 ' 1 ' 

with a stationary temperature 

dT{d/2)i 2 3_ s ^ 2/3 _ 



2tt— (1 - a 2 )ncr d - 1 



1 



16 



It has recently been shown that, in a homogeneous situation, the initial condition is "forgotten" before the stationary 
state has been reached |35|. The system approaches the stationary state through a universal route in which all the 
time dependence of the distribution function, /^(v,*), goes through the instantaneous temperature, T#(i) (in the 
following we will denote the universal state with the subindex H). In contrast with the free cooling case, due to the 
parameter £ 2 , we can construct an additional quantity with dimensions of temperature apart from the instantaneous 
temperature, namely the stationary temperature given by eq. (|13p . Then, by dimensional analysis the distribution 
function has a two parameter scaling form 



Mv,t) 



v H (t)' 



fX(c,/3), 



v H (ty 



v H (t)' 



(14) 



Let us remark that we have redefined the variable c compared to the one introduced in cq. (|10[) and we have 
introduced the instantaneous thermal velocity, vn(t). In reference 35] it was shown that this state actually exists 
and the dynamics is partitioned in a first rapid stage where initial conditions matter, and a subsequent universal 
relaxation towards stationarity, where only the distance to the steady state is relevant, through the dimensionless 
inverse typical velocity /3 = v s /vu{t), i.e. 



/(v,*|/o) 



«&(*) 



x(c,/3) —>/.(«). 



(15) 



Let us note that a similar two parameter scaling occurs in the uniform shear flow of granular gases [371 . |38( . In this 
case, after a quick transient, the system forgets the initial condition and evolves to the stationary state through a 
normal state, where the role of (3 is played by the dimensionless shear rate, a* = a[na d ~ 1 vo(t)]~ 1 . As in the stationary 
state, numerical simulations show that the scaled distribution is close to a Maxwellian and it can be calculated in the 
first Sonine approximation 



x(c,P)~XM(c){l + a 2 (f3)S 2 (c% 



where, by definition, we have 



dcx(c,P) = 1, 



Jdccx(c,l3) =0, J 



dcc 2 X (c,f3) = -. 



(16) 



(17) 



Introducing the approximated distribution (|16[) into the Boltzmann equation, and neglecting the non-linear terms in 
a,2(l3), it is possible to identify the universal distribution, i.e. the universal a 2 {P), to be characterized by [35j 



1 



1-/3' 



B 



1 



l aJ H~3' 1; 



4B 



(18) 
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valid for < < 1 and 

ABa% (1 AB 10 1 \ 

for j3 > 1. Here we have introduced the hyper-geometric function, 2 _Fi, [23| and the coefficient 

73 + 8d{7 - 3a) + 15a[2a(l - a) - 7] 

~~ 16(1 - a) (3 + 2d + 2a 2 ) + af[85 + d(30a - 62) + 3a(10a(l - a) - 39)] ' ^ ' 

As in the long time limit the distribution tends to the stationary state we have 

lira o 2 (/3) = lim a 2 (/3) = a|. (21) 

We also have that the first derivative is continuous 

^1^=^-1^ (22) 

an important property that will turn to be needed in the following sections. 

The evolution equation for the temperature (or equivalently for the thermal velocity, vjj) in the universal state can 
be calculated by inserting the scaling form (|14p into the Boltzmann equation and taking the second velocity moment. 
When this is done, we obtain 

^|W=r 1 [Mi)/3 3 -M«W(i) 2 , (23) 

where we have introduced the dimensionless coefficient 

M0) = ~^ / dc i / rfc 2X (c 1 ,/3)x(c 2 ,/3)T (c 1 ,c 2 )( C 2+ c 2). (24) 

The operator To is 

T (cx,c 2 )= fd&e(c 12 ■ <r)(c 12 ■ a)(b a - 1), (25) 
where b a replaces the velocities vi and v 2 by the postcollisional ones and v 2 given by 

Vi = vi — [a ■ vi 2 )cr, 

. 1 + a 

v 2 = v 2 H — (cr • vi 2 )<t. (26) 

By inserting the function \ hi the nr st Sonine approximation, eq. (|16[) , an approximate evolution equation for wjj can 
be obtained. Nevertheless, we will be only interested on situations where we are close to the stationary state, /3 = 1, 
and then equation (|23|) can be linearized obtaining 

dv H (t) v s 

df = ~7y m(t) - v s ] , (27) 
where we have introduced the dimensionless coefficient 



7 = 3^(1) - 



(28) 

0=1 



d;3 

The coefficient 7 is calculated in Appendix [X] in the first Sonine approximation. Let us note that the difference of eq. 
(j27|) with the equivalent one in [23| is just the term , that, although small, is shown in the Appendix to be 

of the same order than a 2 . The solution of equation (|27|) is 

v H (t)=v s + [v H (0)-v s ]e-^ t , (29) 

that shows that the relaxation of vu to the stationary value is given in terms of the coefficient 7. 

The universal state that we have identified, that in the following will be referred to as the /3-state, is the analogous, 
for heated granular gases, of the homogeneous cooling state for unforced systems 0- It represents the most general 
homogeneous hydrodynamic state and, as will be seen in the next sections, its existence is crucial in the study of the 
relevance of a hydrodynamic description for granular gases heated by the stochastic thermostat. 



5 



III. THE LINEARIZED BOLTZMANN EQUATION 

Let us consider states close to the stationarity, so that we can write 

Sf(x 1 ,t) = f(x 1 ,t)-f s (v 1 ), \Sf( Xl ,t)\ </ s K). (30) 

The evolution equation for this function is obtained from the Boltzmann equation ([3]) by neglecting the non-linear 
terms in Sf, obtaining 

—Sf(x 1 ,t) + v 1 ■ —8f(x u t) = K(yi)6f(xi,t), (31) 
at or i 

where K (vi) is a linear operator defined by 

K(v 1 )=a d - 1 j dv 2 f (v 1 ,v 2 )(l+P 12 )/ s (w 2 ) + ^^2, (32) 

with Pi 2 an operator that interchanges the label 1 and 2 in the function on which it acts. Now, let us introduce a 
dimcnsionlcss space variable as 

l=j, £=(na d - 1 )-\ (33) 
where I is proportional to the mean free path, and a dimensionless time 

s = jt, (34) 

which is proportional to the number of collisions per particle in the interval (0,i). It is also convenient to introduce 
the dimensionless distribution, 8\ through 

TL 

6f(x,t) = - s 5 X (l,c,8). (35) 



The equation for S\ reads 

o r o ■ 

6 X (l,c,s). (36) 



^-Sx{l,c s) = 

OS 



A / \ 9' 



Here we have introduced the homogeneous linearized Boltzmann collision operator, A, defined by 



A(ci)= / dc 2 f (ci,c 2 )(l + Pi 2 )x a (c2) + V«^2. ( 37 ) 



^Q2 

2 del 



~ 9 tit tt— (1 - a 2 ) . 



where 

v/2dr(d/2) V" ^ 16 

is the dimensionless amplitude of the noise calculated in the first Sonine approximation. Since the equation is linear 
and the linearized Boltzmann collision operator does not change the space variable, it is convenient to introduce the 
Fourier component 

<5 X k(c,s) = J dLe- ik l S X (l,c,s). (39) 
For an infinite system or with periodic boundary conditions, the evolution equation for these components is 

£$Xk(c, s) = [A(c) - ik • c] <5 Xk (c, s). (40) 

OS 

Equation (|36|) or the Fourier counterpart, cq. (|40[) , is the so-called linearized Boltzmann equation and it describes 
the dynamics of any small perturbation around the homogeneous stationary state. This provides us with our starting 
point for the study of the possibility of a hydrodynamic description close to the stationary state. 
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The solution of eq. (j40|) can be written formally as 



(41) 



that shows clearly that the excitations of the gas are determined by the spectrum properties of the linear operator 
A(c) — ik ■ c. This suggests the study of the eigenvalue problem 



[A(c)-ik.cfo(k,c) = A 3 -(k)fc(k,c), 
that is posed in a Hilbert space of functions of c with scalar product 

(g\h) = ( dcx- 1 (c).g*(c)%), 



(42) 



(43) 



with g* the complex conjugate of g. The search for hydrodynamic excitations, which are defined as the ones associated 
to the slowest modes, can be carried out by assuming that the modes are analytic in k and looking first for the k = 
solution of eq. (|42|) . So, let us consider the homogeneous eigenvalue problem 



A(c)&(c)=A^(c), 



(44) 



where we have introduced the notation £j(k = 0,c) = £j( c ) and Aj(k = 0) = Xj. We will now see that the special 
universal solution studied in the previous section allows us to identify a family of exact solutions of the linearized 
Boltzmann equation related to d + 2 modes of the homogeneous linearized collision operator, A. The idea is similar 
in spirit to the one introduced in [8l4Iol| to calculate the hydrodynamic eigenfunctions in the free cooling case. Let us 
consider the family of exact solutions of the homogeneous non-linear Boltzmann equation 



f H (v,t) 



v H {t) d 



\ 



VH(t) ' VH(t) 



(45) 



which are parameterized by the density, n, the constant velocity flow, u, and the thermal velocity, vn(t). The bars 
on the quantities refer to variables that differ from (n, vh) introduced earlier, vn{t) being the thermal velocity of the 
/3-state corresponding to n. If we consider states close to stationarity, the function vjj(t) is known and given by eq. 
But then, the family 



*/(v,t) 



v H (t) d ' 



u 



vn{t) ' v H {t) 



fs{v), 



(46) 



to linear order in the fields, has to be a solution of the linearized Boltzmann equation. In Appendix [B] Sf is calculated 
to linear order, which leads to the corresponding scaled distribution 



6x(c,s 

5n 
'■in 



6n 
n 

Sv H (0) 



Xs(c) 



ld_ 
3dc~ 



i c Xs{c)} 



d_ 

dc 



Xs(c) 



13=1 



(47) 



where Sn = n — n and 5vh(0) = vh(0) — v s . The functions given by (|47[) generate a family of solutions of the 
homogeneous linearized Boltzmann equation, that can be seen as the superposition of d + 2 modes. We consequently 
have identified the following eigenfunctions of A 



1 d 

6(c) = Xs(c) + • [cx»(c)\, 
d 

£2( C ) = -^-Xs(c), 

oc 



6(c) 



9c 



i c Xs(c)} 



df3 



x(c,P) 



(48) 
(49) 
(50) 



13=1 



with the corresponding eigenvalues 



Ai = A 2 = 0, 



(51) 
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where the null eigenvalue is (d + l)-fold degenerate. Moreover, cq. (|47|) can be rewritten as 



3 



5 X (c,s) =^e A ^fe(c)| ( 5 X (c,0))^(c), (52) 
i=i 

with the following definitions of the functions £j 

f 1 (c)=Xa(c), l 2 (c)=X,(c)c £ 3 = X.(c) - 0- (53) 

These functions are just the linear combinations of the the first two velocity moments, {1, c, c 2 }, normalized to enforce 

(ti\ti)=8ij, i,j = 1,2,3. (54) 

In fact, the functions £i and £ 2 have been previously identified in (30l . l3ll |40| . where they were used to study the 
fluctuations of quantities like the total energy in the stationary state. There, it was also proven that £i and £ 2 are 
eigenfunctions of the adjoint operator of the linearized collision operator, A + , associated to the null eigenvalue. This 
follows from the conservation of the number of particles and momentum: 

J dcA(c)/i(c) = J dcciA(c)h(c) = 0. (55) 

Then, if we assume that eigenfunctions of A, form a complete set, in the expansion 

oo 

/(c) = E c '^( c )' ( 56 ) 

we have C\ — and C 2 = (£2 1/)- With the aid of eq. (|4"T)) . we are able to identify a new eigenfunction, £3, that, 

as seen in (|50[) . depends on the derivative with respect to /? of the universal state x m the stationary regime. Let 
us remark that £3 is not an eigenfunction of A + . We only have C3 = (£3!/) if the function / belongs to the space 
generated by the set {£,}f =1 . 

In the elastic case, i.e. a = 1 and £q = 0, the spectrum of A(k, c) has been analyzed in detail and it is known 
that the modes associated to the locally conserved quantities, the hydrodynamic modes, are the slowest ones in the 
k — > limit, where the respective eigenvalues vanish [4l|. Furthermore, it is known that the spectrum is analytic in 
k around k = 0, and that there is scale separation, i.e. the hydrodynamic modes are isolated from the rest of modes. 
For the inelastic linearized collision operator, no such result is available. We have nevertheless shown that d+ 1 of the 
d + 2 identified modes are associated to the vanishing eigenvalue, and the remaining mode is associated to 7, which 
itself vanishes in the elastic limit. We therefore expect that a similar property will hold with the identified modes at 
least close to the elastic limit- and these modes will henceforth be coined 'hydrodynamic'. In the following, we will 
assume that they are the slowest ones and that they are analytic. Under this proviso, the asymptotic behavior of the 
one particle distribution function for k <C 1, s ^ 1 is 

d+2 

where the {Kj}^ 2 depend on the initial condition. Note that, as seen in eqs. (|46f and (|47|) . the hydrodynamic 
eigenfunctions given by eqs. (|4"8")) - (j5T)|) , expand the subspace of functions generated by the difference of a local /3-state 
with the stationary state. Then, eq. (|5T[) can be rewritten in the original variables as 



5f r,v,t « x 
v(r,t) a 



v-u(r,i) v s [n(r,t)} 



v(r,t) v(r,t) 



fs(v)- (58) 



This shows that for small gradients (or equivalently for k <C 1) and in the long time limit, i.e. in the time in 
which the non-hydrodynamic modes have decayed, all the time dependence in the distribution function is through the 
hydrodynamic fields. Moreover, the distribution function takes the form of a local /3-state distribution, which plays, 
in this context, the role of a reference state. To evaluate explicitly the time evolution of the fields, it is necessary to 
calculate {Aj(fc)}f+ 2 . If we assume that they are analytic in k 7 this can be carry out studying the eigenvalue problem, 



eq. (|42j) . by standard perturbation theory, as was performed in [8j, |10| for the free cooling case. In the next section, we 
will evaluate the eigenvalues to Navier-Stokes order, i.e. k 2 order, but by a different method, deriving the evolution 
equations for the linear deviations of the hydrodynamic fields. 
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IV. LINEAR HYDRODYNAMICS AROUND THE STATIONARY STATE 



The objective of this section is to derive evolution equations for the deviations of the hydrodynamic fields around its 
stationary values at Navier-Stokes order. The analysis of these equations will clarify the differences and analogies with 
respect to the elastic case. The hydrodynamic eigenvalues, {\j(k)}j=i, will be obtained by identifying the asymptotic 
behavior of the solutions in the proper limit. Let us start writing the deviations of the hydrodynamic fields as 

n(r, t) — n 



p(l,s) = nr ' t — VL = [ dcSx(l,c,s), 
n J 

w(l,s) = U = / dcc6x(l,c, s), 

v s J 



(59) 
(60) 
(61) 



that can be expressed in terms of the scalar products of the distribution function with the functions In 
Fourier space, they read 



(£i(c)|£xk(c,s)) = Pk{s), 
(l2( c )I^Xk(c,s)) = w k (s), 

(6(c)|5 X k(c, S )) = ^ k ( S ) + ^ k ( S ). 
Now, let us define the relevant projector in the hydrodynamic subspace 

d+2 

PMc)=£&(c)|M<0>fc(c), 



and also the orthogonal one 



Q = l-V, 



(62) 
(63) 

(64) 



(65) 



(66) 



where X is the identity operator. As alluded to above, (£,3\h) is not the actual component of h into £3 but, still, eq. 
(|rJ5j) defines a projector as V 2 = 1 (we also have Q 2 =1). 

If we apply the projectors V and Q to the linearized Boltzmann equation, eq. (|40p . we obtain the following set of 
coupled equations 



d_ 

da 

d_ 



- P[A(c) - ik ■ c]V 

- C[A(c)-ik-c]fi 



P<5 X k(c, s) = [?A(c) - ViV ■ c] Q<5 X k(c, s), 
Q<5xk(c, s) = -Qik ■ cVSxk(c, s), 



(67) 
(68) 



where we have used that QAV = 0. However, VAQ ^ because £3 is not a left eigenfunction of A. In fact, we have 



TA(c)Qh(c) = fc(c)<&(c)|A(c)eM<0>. 
The d + 2 components of eq. (|67[) are the evolution equations for the hydrodynamic fields 

8 



ds 
d 



Pk(s) + ik ■ Wk(s) = 0, 



, w k (s) + [p k (s) + k (s)] + ik ■ n k (s) = 0, 
ds 2 



|a(*)+7 
us 



3pk(s)+^k(s) 



2i 2i 

—k ■ w k (s) + — k • k (s) = 5( k (s), 



where we have introduced the pressure tensor and heat flux 



tf k («) = J dc*£(c)Q8 Xk (c,s), k (s) = J dcS(c)Q<5 X k(c,s), 



(69) 

(70) 
(71) 
(72) 

(73) 
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with 



Ajp( c ) — CjCp ~^^jpi ^j( c ) — ( c 



d + 2 



' J ! 



and the deviation of the cooling rate 



6{k(s) = J dc^-A(c)Q5 X k(c, s). 



(74) 



(75) 



Let us note that the operator Q can be skipped in eq. (|73|) because, due to symmetry properties, we have 

(x(c)A J> (c)|^(c)) = (x(c)E,-(c)|^(c))=0 j,p = l,...,d, = 1,... d + 2. (76) 



To evaluate the hydrodynamic equations to k 2 order, we need the fluxes given by eq. (|73|). and consequently Sxu 
to first order in k. This is evaluated in Appendix [Cl for an initial condition of the form QSxk{c,0) = and neglecting 
all the k contributions in the kinetic modes, yielding 



<2<5Xk(c, s) « - / ds'e QA ^ s - s "> Qik ■ cVd X k(c, s'). 
Jo 



(77) 



Note that, as the hydrodynamic eigenfunctions expand the subspace of functions generated by the difference of a local 
/3-state with the stationary state, the condition Q(5xk(c, 0) = represents an initial condition of the local /3-state 
form. When (|77p is inserted in (|73[) . taking into account the symmetries of the system as discussed in Appendix [Cl 
the following expressions for the fluxes are obtained 



nkjp(s) — —i / ds'G xy (s - s') k 3 Wk. p (s') + k p Wk,j{s') - -Sjpk ■ xvu(s') 
Jo L " 



and 



4>k,j(s) = ~ikj J ds' |/9k(s') 
where we have introduced the "correlation" functions 



H 1 (s-s') + -H 3 (s-s') 



+ 2 e *(s')H 3 (s - s') 



and 



G xy (s) = J dcA xy (c)e A{c)s c x & y {c) 
flj-(s) = J dc^ x ( C )e A ^ s c x ^( C ), j = 



(78) 



(79) 



(80) 



1,3. 



It is important to remark that, as A xy and T, x arc orthogonal to the hydrodynamic modes, the functions G xy (s) and 
Hj(s) decay with the kinetic modes. Under the same hypothesis as for the fluxes, the cooling rate, (5£k, is evaluated 
in Appendix [Dl to k 2 order, with the result 

5(k(s) = -2i[ ds'k • w k (s')Z(s - s') 
Jo 

-2k 2 [ ds' (pk(s') 



Z x {s - s') + -Z 3 {s - s') 



+ -e k (s')Zs(s-s')\. 



We have introduced here the functions 



and 



Z{s) = (&(c)|[A(c) - A 3 ]e A ( c ) s c,6,,(c)), 
Zj(s) = (6(c)|[A(c) - Aa] ^ d S 'e A ^ s -^c x Qe A ^ s 'c x ^(c)), j = 1,3, 



(82) 
(83) 

(84) 
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that, due to the bi-orthogonality condition, cq. (|54[) . also decay with the kinetic modes. Note that, at variance with 
the free cooling case Q, there is here a first order in k contribution to the cooling rate. 

It also proves convenient to introduce the parallel and transversal components of the velocity 

w k ,l|(s)=k-w k (s), wgx = #l ■ w k(s), i = l,...,d-l, (85) 

where k = k/fc is a unit vector in the direction of k and {k*f 'j^Zj 1 is an orthogonal basis of the subspace orthogonal 
to k. In terms of these components, the hydrodynamic equations arc 



d_ 

ds 



— pk(s) +ikw kj \\(s) = 0, (86) 



^- w ^(s) + k 2 J^ ds'G xy (s-s')w^ ± (s')=Q, j = l,...,d-l (87) 
§; W K\\( S ) + \ k W{s) +6*(s)] + 2^-^-k 2 j ds'G xy (s-s')w kt \\(s') = 0, 



d 

\(s)+7 



\pki,a) + 6k(a) 



2i 

+ — ku)k,\\ (s) + 2ik / ds'Z(s — s')tUk,||(s') 



-k 2 \ ds'[G 1 (s-s')p k (s') + G 3 (s-s')6 k (s')]=0, (89) 



o 

where we have introduced the new correlation functions 

G 1 (s) = H 1 (s) + ±H 3 (s) + dZ 1 (s) + ^Z 3 (s), G 3 (s) = ±H 3 (s) + ^Z 3 (s). (90) 

It is noteworthy that, apart from the assumption of an initial condition in the hydrodynamic subspace, the only 
approximations made in the derivation of eqs. (|86[) - (|89[) are the expansion to second order in the gradients of the 
hydrodynamic fields and the neglect of the k contribution in the memory kernels that all decay with the kinetic modes. 
Of course, the kernels are for the moment unknown, but we will see later that they can be calculated approximately. 
Before doing so, we evaluate the asymptotic behavior of eqs. (|86|) -()89 j) in the hydrodynamic limit (under which eq. 
([57]) was derived). 

Due to the structure of the equations, it is convenient to introduce the Laplace transforms of the fields 



1 



f(z) = / dae-"f(a), f(s) = — / dze*'f(z), (91) 

JO Z7Tt Jc-ioo 

where c is bigger than the real part of all the poles of /. In the Laplace space the convolutions transform into products 
and Eqs. (|55 ]) -(|g9" ]l become 

zw^(z) + k 2 G xy {z)w^ ± (z) = w^ ± (0), j = l,...,d-l. (92) 

Pk(z) \ ( Pk(0) 
[zI + A(k,z)} | w Kll (z) = ^ k . M (0) | , (93) 

9w(z) I I k(O) 



where we have introduced the matrix 



ik 



A(k,z)=\ ifc 2^G xy (z)k 2 Ik |, (94) 

| 7 +|G 1 (z)fc 2 iq(z)k 1 +l l G 3 {z)k 2 

with 

q{z) = -[l + dZ{z)}. (95) 
d 

The equation for the transverse velocity, eq. (|92[) . is analyzed in detail in Appendix lEl Assuming that G xy (s) is a 
linear combination of kinetic modes (as is expected to be) , we obtain in the hydrodynamic limit 

"ftW^i^', (96) 
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where we have introduced the shear viscosity 

/>oo 

r?= / dsG xy (s). (97) 
Jo 

In the same limit, the solution of the coupled hydrodynamic equations, eq. (|93p . is analyzed in Appendix IFl obtaining 

|y(M) » ^(^|y(k,0)) e ^^|^ 0) ), (98) 
0=1 

where we have introduced the notation 

/ Pk(s) \ 

(99) 



Pk(s) 

|y(k,s)) = I w k ,||0) 



The functions {|^ 0) )}| = _ x are the zeroth order in k contribution to the expansion of the eigenfunctions, \ipp(k, z)), of 
the matrix A(k, z) 

\Mk, *)) = \$\z)) + mf{z)) + k 2 \^;\z)) + ..., 

with A(k, z)\ij)p(k, z)) = ap{k, z)\ijjp(k, z)) . They are calculated in Appendix IFl obtaining 



(100) 



i^ 0) > = I ve i , iv>n = I ve i , ivn = I I , 



-6 



<,<°)\ - 



(0)\ _ 



(101) 



where we see that they do not depend on z. The other set of functions, {(<fip\} 3 R =1 , is the bi-orthogonal set 



il = ( ~.r^,0l , 



12' 2^' 



. 1 1 

y 2| = 777, TT-^ 



12' 2y/6' 



3,0,1), 



(102) 



that is constructed to have (0,3 IV^) = $p,/3'- The scalar product (u\v) between two vectors \u) and \v) is the usual 
Euclidean scalar product. There is no confusion with the one introduced in eq. (|43p because they appear in different 
contexts. The explicit expressions of the hydrodynamic eigenvalues, {\p(k)}^ =1 , to k 2 order are 



Ai(fc) = 




A 2 (fc) = 




As(fc) = 


-7 + 



d 


- 1 








d ?H 


- — 




~d 


- 1 








d H 


- — 

4 7 





k 2 



1 1 2 , , 
3^ + ^'-d (K + Ce) 



(103) 
(104) 
(105) 



(106) 



where the shear viscosity, 77, is given in (|97|) . and we have introduced the heat conductivity, k, as 

1 r 00 

The expressions for the other transport coefficients are 

2 f°° 

qo = - + 2 J dsZ(s), (107) 



q w = 3 + 2 
a 



dse~< s Z{s), 
dse ls Zz{s). 



(108) 
(109) 



Eqs. (|M|) and (j9"5)l with the expressions for the transport coefficients, eqs. (f9"T)) . ()106|) - (|109|) are the main results of 
the paper. It is important to remark that the equations are valid in the linear regime close to the stationary state, 
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but the transport coefficients depend on the structure of the non-stationary /3-state through the eigenvaiue 7 and the 
third cigenfunction. Actually, similar effects are shown to be essential to understand the dynamics of a homogeneous 
perturbation of the temperature close to the steady uniform shear flow for granular gases [421 ■ Let us also note that, 
although there is a coupling between the heat flux and the density as it is seen in eq. (|79p , this coupling is not reflected 
at the level of Navier-Stokes linear hydrodynamic (there is no contribution from H\ (s) in the transport coefficients) . 
This is also the case in the free cooling case where the hydrodynamic eigenvalues to k 2 order do not depend on the 
diffusive conductivity Hoj|. Moreover, the memory kernel Z(s) appears in qo and q w weighted in different ways, 
reflecting the non-Markovian character of eqs. ((55)) - ([5§)l . Nevertheless, as we shall see in the remainder, these effects 
are expected to be small. We also emphasize that the viscosity and heat conductivity have been calculated in [34[ 
applying the Chapmann-Enskog expansion to the inelastic Enskog equation, obtaining equivalent expressions for both 
transport coefficients in the low-density limit. 

Our goal is now to calculate all the correlations functions, G xy {s), G\{s) 1 G 3 (s) and Z(s), that appear in the 
hydrodynamic equations, eqs. (|86|) - (f89|) . in an approximate way. With this, we will be able to obtain explicit formulas 
for all the transport coefficients. The idea is reminiscent of that used for free-cooling systems |43[ and consists in 
treating the functions £3(0), Xs (c)Aj P (c), and \s ( c )Sj(c) as if they were eigenfunctions of the adjoint linearized 
Boltzmann operator, A + (the adjoint is taken with the scalar product of eq. (|4U)) ). That is, we assume 

A + (c)e 3 (c) « A 3 &(c), (110) 

and 

A+(c) Xs ( C )A JP (c) » \$ HXs (c)A jp (c), A+(c) Xs ( c )S,(c) » A^ Xs (c)£,(c), (111) 

where and \ff H are two non-hydrodynamic (kinetic) eigenvalues that have to be calculated consistently. Within 
this approximation we trivially have 

Z(«)«0, Z,-(a)«0, j = 1,3, (112) 

so that 

2 2 

Qo ~ 3, q w ~ j, Ce ~ 0, (113) 
a a 

and 

G 1 (s)^H 1 (s) + ^H 3 (s), G 3 (s)^^H 3 (s). (114) 
Assuming the property ([Hip holds, the functions G xy (s), Hi(s), and H 3 (s) follow by straightforward calculations: 

G xy (s) = ie>&., H l(s) = JA±m±^L ^ (115) 



H 3 (s) = 



/3=1 



(116) 



where the eigenvalues A^^ and A^^ are calculated in Appendix O in the approximation of eq. (|llll) . With these 
functions, the transport coefficient are easily calculated 



1 



2IA 



(i) I 



2 |A}&|-7 



d + 2 da 2 {P) 



d/3 



(117) 



In Fig. [T] the reduced viscosity j s plotted as a function of the inelasticity for d 



2 and d = 3, and in Fig. 



[2] the same is done but for the reduced conductivity _ Although strictly speaking, we do not consider the exact 
same thermostating mechanism as in [26j j where the driving amplitude is chosen to depend on the local temperature, 
it is nevertheless relevant to compare our predictions to those of Ref. [2(1]. In the case of the shear viscosity we 
obtain very similar results. The difference is due to the approximate method to evaluate the coefficient, yielding a 
smoother curve with the present method. In contrast, for the reduced conductivity, we obtain strong discrepancies. 
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0.4 0.6 

a 
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a 



Figure 2. Reduced conductivity of a two-dimensional system (left) and three-dimensional system (right) as a function of the 
inelasticity. The solid line is the theoretical prediction given by eq. 1)1171) . the dashed line is the theoretical prediction of [26l |. 
and the dots are the simulation results of [331 . 



Equation (|117|) predicts an enhanced of the conductivity as the inelasticity increases while the prediction of [26| goes 
in the opposite direction. For completeness, we have reported some simulation data available in the literature. For 
the shear viscosity, the results were taken from Ref. [26[ for d = 3. The agreement with the theoretical prediction 
is good. For the heat conductivity, we have taken the results of [32| . pertaining also to dissipative hard spheres, for 
the smallest density and wave vector. The agreement with equation (JTTTJ) in this case is very good for the two values 
of the coefficient of normal restitution. Let us note that the transport coefficients were also evaluated for moderate 
densities i n 1 441) via the Enskog equation and the density dependence agreed qualitatively well with the simulation 
results of 13211 . 



V. CONCLUSIONS AND PERSPECTIVES 



In this paper we have derived the Navier-Stokes hydrodynamic equations for a system of hard particles heated by 
the so-called stochastic thermostat. We have restricted to situations close to the homogeneous stationary state that 
the system reaches in the long time limit. Under these conditions, the system is described by the Boltzmann equation 
linearized around the stationary state. We could calculate the eigenvalues and eigenfunctions of the operator describing 
the dynamics of the system that are relevant in the hydrodynamic description. Let us remark that, although we are 
considering linear response around the stationary state, the modes depend on the properties of the time-dependent 
/3-state through quantities related to a x(c,ff) 



op 



. The properties of the /3-state brought to the fore were summarized 

3=1 
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in section [TT] and in particular, j3 — 1 can be viewed as measuring the distance to stationarity. With the aid of 
these modes, we derived the linear Navier-Stokes equations obtaining formulas for the transport coefficients that are 
expressed as Green-Kubo relations. Assuming that the time correlation functions that appear in these formulas decay 
with only one kinetic mode, we calculated explicitly the transport coefficients as functions of the inelasticity, a, and 
the spatial dimension, d. Let us note that, at this level, the dynamics also depends on the properties of the /3-state, 
through the transport coefficients. Moreover, as it is reflected in eq. ([58]) . the /3-state plays the role of reference state 
in the sense that, in the hydrodynamic scale a local /3-state distribution is reached. This fact is in connection with 
the results of [l4[ where it is seen that, close to a stationary state, the zeroth order in the gradients distribution is not 
merely the local stationary distribution but a more complex one (in our case played by the /3-state). In this sense, 
for situations where non-linear effects are important, the complete hydrodynamic equations at Navier-Stokes order 
could be derived by the Chapman-Enskog scheme taking into account these ideas. We hope this work will contribute 
to stimulate more studies in this direction and also at the level of computer simulations by, for example, measuring 
the transport coefficients. Also, since many effects shown in the paper depend on the two-parameter distribution of 
the /3-state, it is expected that similar phenomena will occur for other sorts of homogeneous thermostats. 
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Appendix A: Evaluation of 7 in the first Sonine approximation 



In this Appendix we calculate the coefficient 7 given by eq. 

dfi(8) 



7 = 3 M (1) 



dB 



(Al) 



0=1 



in the first Sonine approximation. By substituting the approximate expression of x(c, 8) given by eq. (|16[) into the 
definition of (i(B) (see eq. (f24|) ). we obtain 



where 



and 



fi(B) = fiM + Msa 2 (/3) 



MM = / dCi I dC2XM{ci)XM(c 2 )T (ci,C 2 )(cl + Cj), 



Ms = -~ J dc x J dc2XM{ci)xM(c2)S 2 {c 2 1 )To(c 1 ,C2)(cj + c\). 
The integrals can be readily performed with the result 



MM 



7T— (1 - a 2 ) 
V2dT{d/2) ' 



Ms 



a— 1 , o , 

3tt— (1-q 2 ) 
16y/2dT(d/2) ' 



Finally, the expression for 7 is 



7 = 3flM + 



3o 2 (l) 



da 2 {8) 



dd 



(9=1 



Ms, 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



where the expression of a2(/3) has been provided in the main text. In figure [31 7 is plotted as a function of the 
inelasticity. The approximate expression of the reference [22| . j v n — 3/i(l), is also plotted, finding very similar 

results. Nevertheless, let us note that the difference between the two expressions, 7^^ — 7 = da ^i^ Us is of the 

ap 13=1 

order of 02(1). 
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Figure 3. Eigenvalue 7 as a function of the inelasticity for a two-dimensional system (left) and three-dimensional system (right). 
The solid line is the theoretical prediction for 7 and the dashed lined is the one of [23], j v n. 



Appendix B: Evaluation of the hydrodynamic eigenfunctions 

Let us evaluate the expression 



Sf(v,t) 



v H {ty 



[X 



vn{t) v H (t) 



/.(«), 



(Bl) 



to linear order in Sn = n — n, and Svh(0) = vh(0) — v s . We consider the case u = 0, for assuming u =/= is a 
straightforward generalization. Let us rewrite eq. (|B1[) as 



*/(v,t) 



v H (t) d 



X 



VH(t) ' Vu(t) 



n 




n 








+ *?* 





fs(v). 



(B2) 



Then, the first two terms on the left hand side of the equation are the difference between the /3-state and its corre- 
sponding stationary state (both are linked to the same density, n), while the last two terms are just the difference 
between two very close stationary states. For the last two terms and, taking into account that 



as follows from from eq. (|13p which defines the unspecified function /C, and thus 



(B3) 



—y f-V-Y (-) =— x {- 



V 



V" 



IC5v s 



(d + 3) 

„,d+4 



K 



Xs(c) 



_L d_ 

dv s Xs U 



(B4) 



where 5v s — v s — v s . It is more convenient to write this expression in terms of the difference in densities. Due to eq. 
(IB3D. we have 



Sn Sv* 
— = -3— 

n v s 



so that 



r,<l 









n 




(i) 








n 1 



Xs(c) + ~ • [c%,(c)] 



Now, let us evaluate the other difference 

v v. 



v H (ty 



rX 



VH(t) VH{t) 



vf Xs \v s J dv H {0) \v H (t) d 



d 



vu{t) ' v H (t) 



[vh(0)-v s ], 



(B5) 



(B6) 



(B7) 
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where the subscript s refers to the functional in the stationary state, i.e. F[n,vu( 
derivative is 



d 



dv H (0) \ v H (t) d 



FX 



y H {t) ' v H (t) 

while the partial derivative can be calculated taken into account eq. 

'dv H (t) 



VH(t)d+ l la --[cx(c/3)]+/3 ¥ x( C ,/3) J q Vh{q) 



F[n,v s ]- The functional 
dv H (t) 



giving 



_dv H (0)_ 

Then, by substituting eq. (|B8|) into eq. (|B7|) and defining c ~ v/v s we have 



v H {ty 



VH{t) ' V H (t) 



-7-f* 



5v H (0) ISn 
Vf, 3 n 



where we have substituted n by n and i> s by v s (that can be done to linear order) and we have used that 

v H (0) -v s +v s -v s _ Sv H (0) ISn 
v s v s 3 n 

Taking into account (|B6j) and ()B10j) . we finally arrive at the result reported in the main text for u = 

Sn 
n 



<^x(c,s) 



1 d 

Xs(c) + • [cxs(c)\ 



Sn 5vh(0) 
3n v s 



(B8) 



(B9) 



(BIO) 



(Bll) 



(B12) 



Appendix C: Evaluation of the fluxes to first order in k 

In this appendix, we evaluate the fluxes given by eq. (|73[) to first order in k. The function QSxk fulfills eq. ([68 
and can be integrated formally as 



Q5Xk(c,s) = e s ( A - ikc ) s Q5 X k(c,0)- f d S 'e QiA ~ lkc ^ s - s '^ Qik ■ cVS X k(c, S '). 

Jo 



(CI) 



Choosing the initial condition in the hydrodynamic subspace, i.e. Q<!>Xk(c, 0) = 0, and neglecting the k contribution 
in the kinetic modes, we have 



Q<5xk(c,s)«- f ds'e QAQ ( s - s '*>Qik-cV6 X k{c,s'). 
Jo 

Inserting the above equation into the expression of the pressure tensor and taking into account that 

J dcA TP (c)QA(c)Qh(c) = J dcA jp (c)A(c)h(c), 



we obtain that 



n kjP (s) « - J dcA jp (c) J ds'e A ( s - s '\k • cVS X *(c, s') 

- / dcA jp (c) / ds'^'-^ik-c^&{c)\6xkM)&{c) 
J Jo 0=1 



(C2) 



(C3) 



(CM) 



r q 
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where we have taken into account that there is no coupling with the density nor with the temperature. Finally, by 
symmetry considerations, we arrive at the expression used in the main text 

f s 2 
IIk,jp(s) w —i I ds'G xy (s — s')[kjV^c lP (s') + k p wu,j(s') - -5j p k-wic(s')], (C5) 



where 



G xy (s) = J dcA xy (c)e As c x & y (c). (C6) 
To evaluate the heat flux to first order in k, we substitute expression (|C2j) into the heat flux and taking into account 

J dcE J -(c)QA(c)Q/i(c) = J dcS 3 (c)A(c)/i(c), (C7) 



we obtain 



Ki( s )~-J dcZjic) ds'e A(W) ik- cPSxk(c,s') 

= - dcEjtc) / ds'e A ^- s 'hk-cJ](^(c)\Sxk(c, S '))^(c) 
J Jo p=1 

= -ikj [ dcE^c) / S d S 'e A ( s - s ') C ,[(ei(c)|^k(c, S ')>ei(c) + (6(c)|<5xk(c, S ')>6(c)] 



(C8) 



where we have used that there is no coupling with the flow velocity. Finally, taking into account eqs. ([62]) and ([64] 
and symmetry considerations, we arrive at the expression of the main text 



0k j O) = -ikj I ds' <{ pk(s') 

where 



Fi(s-s') + iiT 3 (s-s') 



+ le k (s')H 3 (s-s') 



Hj(s) = J dc^ x {c)e As c x i {c), j = 1,3. (CIO) 

Appendix D: Evaluation of the cooling rate to second order in k 

We work out here the cooling rate given by eq. (|75|) to second order in k. As in Appendix [Cj choosing an initial 
condition with Q6xk(c, 0) = 0, and neglecting the k contribution in the kinetic modes, we have to k 2 order in the 
hydrodynamic fields 

QS X k(c, s) « - / ds'e QAQ ( s - s '^ Qik ■ cV6 X k{c, a') 
Jo 

- ds' ds"Qe A{s - s '- s " ) k-cQe As "k-cP8xu(c,s'), (Dl) 

Jo Jo 

where the identity for any two operators A and B 

e (A+B)t = e At + f* dt ' e A(t-t') Be (A+B)t'^ (D2) 



has been used in order to perform the expansion. The first term of the right hand side of (|D1|) is the first order in 
k contribution, QSx\t(c, s)W , calculated in Appendix [C| while the second term is the k 2 contribution, QSxk(c, s)( 2 \ 
By substituting ([Dip into ([75]) the cooling rate to k 2 order is obtained. The first order contribution is 



<5Ck 1} (s) = - J dc^-A(c)^ d S 'Qe Ms ~ s '\k-cV5xk(c, S ') 



-2<&(c)|A(c) I" ds'Qe^-^ik-cPSx^s')), (D3) 
Jo 
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where we have used cq. (1531) . Taking into account that 

(6(c)|A(c)Q/i(c)) = (&(c)|[A(c) - \ 3 }h(c)), 

it can be rewritten as 

5(£\s) = -2(6(c)|[A(c) - A3] f ds'e^ s - s \\L- cW Xk (c, s')) 

Jo 

= -2i Vfc p [ S d S 'w Kp ( s ')(^(c)\[A(c) - A 3 ] e A(s - s ')c P 6, P (c)), 
Jo 



(D4) 



(D5) 



where we have used that, by symmetry considerations, there is only coupling with the flow velocity. Finally, we can 
transform (ID5I) to write it as in the main text 



<5Ck(s) (1) = -2i f ds'k ■ w k (s')Z(s - s'), 
Jo 



where 

Z(s) = <f 3 (c)|[A(c) - \ 3 }e As c x & , x (c)). 
By similar manipulations, the second order in k contribution can be written as 

<5Cl 2) (s) = -2(£ 3 (c)|[A(c)-A 3 ] fds' f S d s " e A ( s -'-")k-cQ e As "k-cW X k(c, s ')) 



(D6) 



(D7) 



Jo 

s—s' 



-2<£ 3 (c)|[A(c)-A 3 ] / ds' I d S "e A ( s - s '- s ")k.cQe As "k.c[ Pk ( S ')6(c) + (6(c)|5xk(c, S ')>e3(c)]> 
Jo Jo 

-2fc 2 (£ 3 (c)|[A(c)-A 3 ] fds' f ' ds'' e A ( s - s '- s '') Cx Q e As '' Cx [p k ( S ')6(c) + (6(c)|5xk(c, S '))C3(c)]>. 
Jo Jo 



that can be recast as 



where 



6£>(s) = -2k 2 / ds'ip^Z^s - s') + (e 3 (c)|*Xk(c, s'))Z 3 (s - s% 
Jo 



Zj(s) = (f 3 (c)|[A(c) - A 3 ] / d S 'e A ( s - s ') Cx Qe As 'c x ^(c)), j = 1,3. 

Jo 



or, as written in the main text 



S(i 2 \s) = -2k 2 I ds'\p k (s') 



Zi(s - s') + -Z 3 (s - s') 



+ -6*{s')Z 3 {s - s') 



(D8) 
(D9) 

(D10) 
(Dll) 



Appendix E: Analysis of the equation for the transversal velocity 



In Laplace space, the transversal velocity is given by eq. 

«'k,i(0) 



i"k,i(z) 



z + k 2 G xy {z) ' 



(El) 



where we have skipped the superscript j. This formula depends on G xy (z) that, in principle, can be a complicated 
object. We already know that in real time it decays with the kinetic modes and in the free cooling case it has been 
shown numerically that it is very well fitted by a single exponential (single mode approximation). Let us consider 
this simple form for the general case can be performed in a similar way. In this approximation, we have 



G xy (s) — Ce s , G xy (z) 



C 
z + X' 



(E2) 
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and cq. (|E1[) is just 



that can be inverted exactly in terms of two exponentials [451 . The function given by cq. (|E3|) has the two poles 



Wk ' lW= ,(z + A)+CF (E3) 



... -A + VA 2 - 4Cfc 2 C 2 
«i(fc) = g ~~A ' ( - 1 



and it can be written as 



Then, we have 



... -A - VA 2 - 4Cfc 2 . , C ;2 
«2(*) = ^ ~ -A + yfc 2 , (E5) 



tSk,± (*) = 4n + 4^ • ( E6 

z — zi(fcj z — Z2{k) 



Wk,±(0)(« + A) = A(fc)[* - z 2 (fc)] + S(fc)[z - 

« A(jfe) f z + A - ^kA + B(k) (z + ^-kA , (E7) 



A / w V A 
with which we identify the constant A and B to zeroth order in k 

Anwu,j_(0), B*0. (E8) 
With eqs. (|E6[) and (|E8[) we obtain the expression for the transversal velocity of the main text, eq. 

™k.i(^%,i(0)e-" t2s , (E9) 

with 

^ poo 

V=t= dsG xy (s). (E10) 
A Jo 

If the function G xy {s) is a linear combination of kinetic modes, the analysis can be performed following the same 
lines obtaining eq. (|E9[) with r\ = J Q dsG xy {s). 



Appendix F: Analysis of the coupled hydrodynamic equations 

We evaluate here the asymptotic behavior of eqs. 

/ p k (z) \ / p k (0) \ 
[zI + A{k,z)] w k ,\\(z) = ufc lN (0) , (Fl) 
V *k(*) / V 0k(O) / 

in the hydrodynamic limit. The matrix A(k, z) is given in eq. (|94[) and can be written as 

A(k,z) = A + tkA 1 (z) + k 2 A 2 (z), (F2) 

where 

\ / 1 

A = I 0, = | \ I , (F3) 



7 7 / V g(z) 





2^±G„,(z) 
|G!(z) \G z {z) 



(F4) 
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In terms of the eigenvalues and eigcnfunctions of the matrix A{k, z) 

A(k,z)\ipp(k,z)) = ap{k,z)\ipp(k,z)), 

the solution can be written explicitly as 



(F5) 



\y(k,z)) = [zI + A(k,z)}- 1 \y(k,0))=J2 



—J z + ap(k,z) 



(F6) 



where we have introduced the notation 



|y(k,s)) = \ w u ,\\(s) 



\y{k,z)) = ( Wk,||0) | , 



(F7) 



and the functions {(ipp(k, z) | }| =1 are the left eigenfunctions of A(k,z). Let us introduce the expansion in powers of 
k of the eigenvalues and eigenfunctions 



\MK z)) = \1>f\z)) + k\1>£\z)) + k 2 \4 2) (z)) +..., 
ap(k,z) = af(z) + kaf(z) + k 2 af{z) + .... 
In the hydrodynamic limit, and applying the same ideas as in Appendix [EJ we obtain 



j(k,z))*Yl 



(^|y(k,0)) (0) 
z-\p(k) Wf> h 



(F8) 
(F9) 

(F10) 



where {A l g(fc)}|_ 1 are the hydrodynamic eigenvalues that appear as the smallest root of eq. z + ap(k, z) = 0. The set 
{(4>p\}p = i is the bi-orthogonal set constructed to have {<t>p\il>pi) = $p,P'- The problem is then to calculate the sets 
{iV'fl )}|=u {(^|}|=i an d {ap(k, z)}p =1 . As {A l s(fe)}| =1 is needed to k 2 order, {ap(k, z)}| =1 have to be calculated 



at the same order, which can be done by standard perturbation theory. 
The eigenvalues of Aq can be easily calculated, obtaining 



(o) _ Jo) 



AO) 



so that the vanishing eigenvalue is two-fold degenerate. The corresponding eigenfunctions are 



(Fll) 



lux) = 



M = 1 



l«3> = l^ 0) > 



(F12) 



and the bi-orthogonal set is 



( Vl \ = ( -,-1,0 



(v 2 \ = (0,1,0), (v 3 \ = (<f> 3 \ = [jj, 



In the non-degenerate case, the first contributions to the expansion are 



and 



4 1} (z) = («3|iAi|u3> =0, 



a (2 \z) = (v 3 \A 2 \u 3 ) + - 5^(«„|iAi|«3>(«a|*Ai| 

^ n^3 



so that we have 



a 3 (k,z) ps 7 + 



(F13) 



(F14) 



(F15) 



(F16) 
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M 




1) 









For the degenerate case, we first have to consider the sub-matrix 

Its eigenvalues are the first order correction of the degenerate eigenvalues 

ai\z) = --j=, 4 1) (z) = -^=, 
and its eigenfunctions give us the components of the corresponding zeroth order eigenfunction 

<,(°)\ \ / 1 n\ ( u/,(o)\ ' 



-1 



1 



(Vl|^ 2 

4' 



1 



1 



so that 



(0)\ 



-6 
V6 



(0)\ _ 
2 / _ 

4 / \ -4 



1 4 



6 

a/6 I , 



that do not depend on z. The corresponding bi-orthogonal functions can also be calculated obtaining 



1 1 



;,0 , 



1 1 



,0 



12'2V6'7' Vl2'2V6 
With these functions, the k 2 corrections to the degenerate eigenvalues follow from straightforward calculations 

"2 



and 



a\ 2 \z) = (&|A 2 |<>) - ^U^f^^A^f) = ^-G xy (z) + ±- 



4-, 



9(2) 



3 



so that we have 



a 2 (k,z) w + 



fc 2 , 
A- 2 . 



(F17) 
(F18) 

(F19) 
(F20) 

(F21) 

(F22) 
(F23) 

(F24) 
(F25) 



The smallest root of z + ap(k,z) = 0, with ap(k,z) given by eqs. (|F24|) . (|F25|) and (|F16|) are the hydrodynamic 
eigenvalues given in the main text, i.e. eqs. (|103[) - (|105[) . 



Appendix G: Evaluation of the kinetic eigenvalues 

Let us assume 

A+(c) Xs (c)A ip (c) « A^ Xs (c)A, p (c). 
Multiplying by c x c y and integrating we obtain 

ASJi, tell- /^c > A + (c)x.(c)cc= f.icocAWx.fOcc, 



so that we have 



- 



h = dcc x c y A(c)xs(c)c x Cy. 



(Gl) 



(G2) 



(G3) 



The heating docs not contribute to the integral 1\ and we have 

h = J dc\ J dC2Cl x Clyfo{ci 1 C 2 ){l + Pl2)Xs{ci)Xs{c2)Ci x Cly. 

Taking into account 

a i\/ \ (1 + a) 2 ,. %2- ~ l + a „ . 

(Oct - IJlcixCiy + c 2x c 2y ) = [a ■ c 12 ) a x (jy — (er • c 12 ){c 12y a x + c 12x a y ), 

and the solid angle integrals 

d—l 

/ da<d(a- ■ ci 2 )(<7 • c 12 ) 3 a x a y = fd +5\ c r2Ci 2x c 12y , 

d-l 

2 7T 2 

d&Q(& ■ c 12 )(a ■ ci 2 ) <j = - , d+3 ^ ci2Ci2j, 



we finally obtain 

h 



3(1 + a) 2 1+, 



7T 2 dC! dc 2 Xs(c 1 )Xs(c 2 )cl x ClyCl 2 Cl 2x Cl 2 y 



d—l 



(2d + 3-3a)(l + a)7r— / 23 e 

1 + ~ a 2 ) ) 



2^2d(d + 2)r(d/2) V 16 

where the last integral has been performed with Xs( c ) m the first Sonine approximation. 

(2) 

To calculate X NH we multiply eq. 

A+(c) Xs (c)Si(c) « xf H X.(c)H j (c), 
by , and proceed with an integration over c to obtain 



X (2) 
N H 



(d + 2)u, 2 

The heating contribution to I 2 is simply 



—h, h = / dcc x c 2 A(c)x s (c)c x 

a 2 J 



where £ 2 is given in (|38|) and the collisional term is 

I 2 C) = J dc x c lx c\ J dc 2 T (c 1 ,c 2 )(l + Pi 2 )Xs{ci)Xs{c 2 )ci a 



Taking into account 



(b„ - l)(c\c lx + c\c 2x ) = ( 1 + a ) (<r ■ ci 2 ) 2 [ci x + C2X + 2(<r • c{)a x + 2(<x ■ 02)63] 



2 

and the solid angle integrals 



2 

(<r • c 12 )[cfa x - + 2(6- ■ ci)cia, - 2(<r • c 2 )c 2x ], 



l + a,- u2- 2 



y d&0{& ■ c 12 ){& ■ C12) 3 = pTgjgj cL 
y d&<d(& ■ c 12 )(& ■ c 12 y (<r ■ ci)a x = 2r / rf+5 x [ci2 c i^ + 3ci 2 Ci 2 x(ci • Ci 2 )J, 
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where 



we have 

(q) (1 + aW - ~ f f 

2 2T( d + T j~ J dCl J dc2Xs ^ Xs ^ ClxF ^ Cl ' C2 ^ ( G16 ) 

1 + a 

•7 7 (ci,c 2 ) = - [c 12 (ci a + c 2;E ) + 3ci 2 ci 2a: (ci + c 2 ) • c 12 ] - (cf - cfje^c^x 
a + 3 

1 + a 

H ^ — C 12( c l^ + c 2^) + 2ci 2 [(ci2 • c 2 )c 2x - (ci 2 • Ci)cia;]. (G17) 

Evaluating the integral in the first Sonine approximation we get 

. d-l 

4 C) = - 32^^/2) 1(32 + 16rf)(1 " a) + 4[7 ° + 4?d ~ 3(34 + 5d)a]} - (G18) 
The eigenvalue can be finally written as 

^ = ra(^ + ' f> > Pi.) 

Let us note that the numerator goes as a 2 in the elastic limit, so that a finite result is obtained in that limit. 
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